Source code for hysop.numerics.interpolation.polynomial

# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.


import itertools as it
import numpy as np
import scipy as sp
import sympy as sm
import warnings

try:
    import flint

    has_flint = True
except ImportError:
    import warnings
    from hysop.tools.warning import HysopPerformanceWarning

    msg = "Failed to import python-flint module, falling back to slow sympy solver."
    warnings.warn(msg, HysopPerformanceWarning)

    flint = None
    has_flint = False

from hysop.tools.henum import EnumFactory
from hysop.tools.htypes import check_instance, InstanceOf, to_tuple
from hysop.tools.cache import update_cache, load_data_from_cache
from hysop.tools.io_utils import IO
from hysop.tools.sympy_utils import tensor_xreplace, tensor_symbol
from hysop.tools.decorators import debug
from hysop.tools.warning import HysopCacheWarning
from hysop.numerics.stencil.stencil_generator import CenteredStencilGenerator, MPQ
from hysop.numerics.interpolation.interpolation import MultiScaleInterpolation


def _check_matrices(*x):
    return all(_check_matrix(xi) for xi in x)


def _check_matrix(x):
    return np.isfinite(np.asarray(x).astype(np.float64)).all()


PolynomialInterpolation = EnumFactory.create(
    "PolynomialInterpolation",
    [
        "LINEAR",  # requires 0 ghosts (no derivatives required)
        "CUBIC",  # derivatives order is specified with SpaceDiscretization
        "QUINTIC",  # derivatives order is specified with SpaceDiscretization
        "SEPTIC",  # derivatives order is specified with SpaceDiscretization
        "NONIC",  # derivatives order is specified with SpaceDiscretization
        "CUBIC_FDC2",  # requires 1 ghosts (estimate derivatives with 2nd order centered fd)
        "CUBIC_FDC4",  # requires 2 ghosts (estimate derivatives with 4th order centered fd)
        "CUBIC_FDC6",  # requires 3 ghosts (estimate derivatives with 6th order centered fd)
        "QUINTIC_FDC2",  # requires 1 ghosts
        "QUINTIC_FDC4",  # requires 2 ghosts
        "QUINTIC_FDC6",  # requires 3 ghosts
        "SEPTIC_FDC2",  # requires 2 ghosts
        "SEPTIC_FDC4",  # requires 3 ghosts
        "SEPTIC_FDC6",  # requires 4 ghosts
        "NONIC_FDC2",  # requires 2 ghosts
        "NONIC_FDC4",  # requires 3 ghosts
        "NONIC_FDC6",  # requires 4 ghosts
    ],
)


[docs] class PolynomialInterpolator:
[docs] @classmethod def build_interpolator(cls, pi, dim, fd=None, verbose=False, approximative=False): check_instance(pi, PolynomialInterpolation) kwds = {"dim": dim, "verbose": verbose, "approximative": approximative} spi = str(pi) if spi.startswith("LINEAR"): deg = 1 fd = 2 elif spi.startswith("CUBIC"): deg = 3 elif spi.startswith("QUINTIC"): deg = 5 elif spi.startswith("SEPTIC"): deg = 7 elif spi.startswith("NONIC"): deg = 9 else: msg = f"Unknown PolynomialInterpolation value {pi}." raise NotImplementedError(msg) for i in range(1, 5): if spi.endswith(str(2 * i)): fd = 2 * i break if fd is None: msg = "Could not determine finite differences order." raise RuntimeError(msg) obj = cls(deg=deg, fd=fd, **kwds) obj.spi = spi return obj
[docs] @classmethod def cache_file(cls): _cache_dir = IO.cache_path() + "/numerics" _cache_file = _cache_dir + "/polynomial.pklz" return _cache_file
def __init__(self, dim, deg, fd, approximative=False, verbose=False): """ Create a PolynomialInterpolator. Parameters ---------- dim: int Number of dimensions to interpolate. deg: int or tuple of ints Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...) Degree should be odd on each axis. fd: int or tuple of ints Order of centered finite differences used to compute derivatives in each direction. Will affect the number of ghosts of the method. Should be even because this interpolator only use centered dinite differences. approximative: bool Use np.float64 instead of exact fractions to compute weights. verbose: bool Enable or disabled verbosity, default to False. Attributes ---------- dim: int Number of dimensions to interpolate. deg: tuple of ints Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...) Degree should be odd: deg=2k+1 fd: tuple of ints Order of centered finite differences stencils used to compute derivatives for each direction. p: tuple of ints Corresponds to deg+1. The total number of polynomial coefficients corresponds to P=p0*p1*...*p(dim-1). P: int The total number of polynomial coefficients P=p0*p1*...*p(dim-1) k: tuple of ints Max derivative order required to compute the polynomial interpolator coefficients in each direction. Also the regularity of the resulting interpolant. Corresponds to (deg-1)/2. ghosts: tuple of ints Return the number of ghosts required by the interpolator on each axis. Corresponds to (k>0)*[fd//2 - 1 + (k+1)//2] deg k (k+1)/2 | FDC2 FDC4 FDC6 linear: 1 0 0 | 0 0 0 cubic: 3 1 1 | 1 2 3 quintic: 5 2 1 | 1 2 3 septic: 7 3 2 | 2 3 4 nonic: 9 4 2 | 2 3 4 n: tuple of ints Corresponds to 2*(ghosts+1), the number of required nodes to generate the polynomial coefficients (in each direction). In total we have N=n0*n1*...*n(dim-1) input nodes. G1 G1 <-> <-> X X X X X P P X X P P X X X X X <-----> n1 N: int Total number of input nodes N=n0*n1*...*n(dim-1). M: np.ndarray Grid values to polynomial coefficient matrix: M.dot(F.ravel()) will give C.ravel(), coefficients of P(x0,x1,...) N <---------> X X X X X X ^ |f0| ^ |c0| ^ X X X X X X | |f1| | |c1| | M = X X X X X X | P F = |f2| | N C = M*F = |c2| | P X X X X X X | |f3| | |c3| | X X X X X X v |f4| | |c4| v |f5| v If approximative is set to True, M will contain np.float64 Else is will contain rationals. See Also -------- :class:`PolynomialSubgridInterpolator`: Precompute weights for fixed subgrid interpolation. """ assert dim > 0, "dim<=0" deg = to_tuple(deg) if len(deg) == 1: deg *= dim check_instance(deg, tuple, values=int, size=dim) fd = to_tuple(fd) if len(fd) == 1: fd *= dim check_instance(fd, tuple, values=int, size=dim) p = tuple(degi + 1 for degi in deg) k = tuple((degi - 1) // 2 for degi in deg) ghosts = () n = () for fdi, ki in zip(fd, k): if ki > 0: gi = (fdi // 2) - 1 + (ki + 1) // 2 else: gi = 0 ni = 2 * (gi + 1) ghosts += (gi,) n += (ni,) check_instance(deg, tuple, values=int, size=dim) check_instance(fd, tuple, values=int, size=dim) check_instance(p, tuple, values=int, size=dim) check_instance(k, tuple, values=int, size=dim) check_instance(n, tuple, values=int, size=dim) check_instance(ghosts, tuple, values=int, size=dim) assert all(fdi % 2 == 0 for fdi in fd), "fd % 2 != 0" assert all(degi % 2 == 1 for degi in deg), "deg % 2 != 1" assert all(pi % 2 == 0 for pi in p), "p % 2 != 0" assert all(ni % 2 == 0 for ni in n), "n % 2 != 0" P = np.prod(p, dtype=np.int32) N = np.prod(n, dtype=np.int32) self.dim = dim self.deg = deg self.p = p self.P = P self.k = k self.n = n self.N = N self.fd = fd self.ghosts = ghosts self.approximative = approximative self.verbose = verbose self.key = ("PolynomialInterpolator", dim, deg, fd, approximative) self._build_interpolator() def _collect_stencils(self): dim = self.dim ghosts = self.ghosts verbose = self.verbose fd = self.fd approximative = self.approximative k = self.k n = self.n ghosts = self.ghosts if verbose: print("\nCollecting 1D stencils:") SG = CenteredStencilGenerator() SG.configure(dim=1, dtype=np.float64) S = {} for direction in range(dim): if verbose: print(f" Direction {direction}") Sd = S.setdefault(direction, []) nd = n[direction] kd = k[direction] fdd = fd[direction] gd = ghosts[direction] for i in range(kd + 1): msg = "Failed to compute stencil derivative={}, order={}, origin={}" msg = msg.format(i, fdd, gd) try: if approximative: Si = SG.generate_approximative_stencil(order=fdd, derivative=i) else: Si = SG.generate_exact_stencil(order=fdd, derivative=i) Si.replace_symbols({Si.dx: 1}) except: print(msg) raise msg += f" got {Si.coeffs}." assert not Si.is_symbolic(), msg Si = Si.reshape((nd - 1,)) assert Si.origin == gd Si = Si.coeffs Sd.append(Si) if verbose: print(f" {i}-th derivative: {Si}") return S def _build_stencil(self, dvec): dvec = np.asarray(dvec) k = self.k S = self.S assert dvec.size == self.dim, "dvec.size != dim" assert all(dvec >= 0), f"dvec < 0 => {dvec}" assert all(di <= ki for (di, ki) in zip(dvec, k)), f"dvec > dmax => {dvec}" Sd = S[0][dvec[0]].copy() for d, i in enumerate(dvec[1:], 1): Sd = np.tensordot(Sd, S[d][i], axes=0) return Sd def _build_interpolator(self): dim = self.dim deg = self.deg k = self.k n = self.n p = self.p ghosts = self.ghosts approximative = self.approximative verbose = self.verbose xvals, xvars = tensor_symbol("x", shape=(dim,)) fvals, fvars = tensor_symbol("F", n, ghosts) pvals, pvars = tensor_symbol("C", p) self.xvals, self.xvars = xvals, xvars self.fvals, self.fvars = fvals, fvars self.pvals, self.pvars = pvals, pvars try: data = load_data_from_cache(self.cache_file(), self.key) if data is not None: (P0, S, M) = data if _check_matrix(M): self.P0 = P0 self.S = S self.M = M return except Exception as e: msg = f"Failed to load data from cache because:\n{e}" warnings.warn(msg, HysopCacheWarning) P0 = 0 for idx in it.product(*tuple(range(0, pi) for pi in p)): P0 += pvals[idx] * np.prod(np.power(xvals, idx)) self.P0 = P0 S = self._collect_stencils() self.S = S if verbose: print("\nGenerating variables:") print(" *space vars: ") print(xvals) print(" *grid values:") print(fvals) print(" *polynomial coefficients:") print(pvals) print(" *polynomial patch:") print(P0) print("\nBuilding system...") eqs = [] for dvec in it.product(*tuple(range(0, ki + 1) for ki in k)): if verbose: print(f" => derivative {dvec}") dP0 = P0 for i, deg in enumerate(dvec): dP0 = sm.diff(dP0, xvals[i], deg) stencil = self._build_stencil(dvec) if verbose: print(" stencil:") print(stencil) for idx in it.product(*tuple(range(gi, gi + 2) for gi in ghosts)): if verbose: print(f" -> point {idx}") pos = np.asarray(idx) - ghosts pos = dict(zip(xvals, pos)) eq = dP0.xreplace(pos) for offset in it.product(*tuple(range(-gi, gi + 1) for gi in ghosts)): fidx = tuple(np.add(idx, offset)) sidx = tuple(np.add(offset, ghosts)) eq -= fvals[fidx] * stencil[sidx] eqs.append(eq) if verbose: print(f" {eq}") # Build system such that A*c = B*f where c are the polynomial coefficients and # f the node values dtype = np.float64 if approximative else object A = np.empty((self.P, self.P), dtype=dtype) B = np.empty((self.P, self.N), dtype=dtype) assert len(eqs) == self.P for i, eq in enumerate(eqs): for j, ci in enumerate(pvars): A[i, j] = +eq.coeff(ci) for j, fi in enumerate(fvars): B[i, j] = -eq.coeff(fi) # C = Ainv*B*f = M*f if verbose: print("\nSolving system...") if approximative: Ainv = np.linalg.inv(A) elif has_flint: coeffs = list(flint.fmpq(x.p, x.q) for x in A.ravel()) Afmpq = flint.fmpq_mat(*(A.shape + (coeffs,))) Afmpq_inv = Afmpq.inv() coeffs = list( sm.Rational(sm.Integer(x.p), sm.Integer(x.q)) for x in Afmpq_inv.entries() ) Ainv = np.asarray(coeffs).reshape(A.shape) else: # /!\ sympy is really slow Ainv = np.asarray(sm.Matrix(A).inv()) if verbose: print("\nBuilding matrix...") M = Ainv.dot(B) self.M = M update_cache(self.cache_file(), self.key, (P0, S, M))
[docs] def interpolate(self, fvals): """Return the polynomial interpolating input node values""" fvals = np.asarray(fvals) assert fvals.shape == self.fvals.shape pvals = self.M.dot(fvals.ravel()) P0 = self.P0.xreplace(dict(zip(self.pvars, pvals))) return sm.utilities.lambdify(self.xvars, P0)
[docs] def generate_subgrid_interpolator(self, grid_ratio, dtype=None): return PolynomialSubgridInterpolator( interpolator=self, grid_ratio=grid_ratio, dtype=dtype )
def __hash__(self): objs = (self.dim, self.deg, self.fd, self.approximative) return hash(objs)
[docs] class PolynomialSubgridInterpolator: def __init__(self, interpolator, grid_ratio, dtype=None): """ Create a PolynomialSubgridInterpolator from a PolynomialInterpolator and a number of subrid points. Parameters ---------- interpolator: PolynomialInterpolator Interpolant used to compute weights. grid_ratio: tuple of int Tuple of integers representing the ratio between the coarse and the fine grid. dtype: np.dtype Force to cast dtype for all matrices (interpolator.M may contain rationals). Attributes ---------- dim: int Number of dimensions to interpolate (same as interpolator.dim). ghosts: tuple of int Number of required ghosts. n: tuple of int Corresponds to 2*(ghosts+1), the number of required nodes to generate the polynomial coefficients (same as interpolator.n). N: int Total number of input nodes N including ghosts (same as interpolator.N). N = n0*n1*...*n[dim-1] s: tuple of int Corresponds to grid_ratio + 1, number of points of the subgrid in each directions. Example for a grid ratio=(3,3), we have s=(4,4): O=coarse grid nodes, X=fine grid nodes Coarse grid: Fine grid: ^ O-----O ^ O X X O | | | | X X X X 1 | | | 4 | X X X X v O-----O v O X X O <-----> <-----> 1 4 S: int Represents the number of fine grid points contained in a coarse grid cell. S = s0*s1*...*s[dim-1] gr: tuple of int Corresponds to grid_ratio, number of points of the subgrid in each directions, minus one. Example for a grid ratio=(3,3), we have gr=(3,3) and s=(4,4): O=coarse grid nodes, X=fine grid nodes, -=excluded find grid nodes Coarse grid: Fine grid: ^ O-----O ^ O X X O ^ | | | gr0 | X X X - | s0 1 | | | v X X X - | v O-----O O - - O v <-----> <---> 1 gr1 GR: int Represents the number of fine grid points contained in a coarse grid cell exluding right most points. GR = gr0*gr1*...*gr[dim-1] W: np.ndarray Pre computed weights to interpolate directly from coarse to fine grid. Let F be the vector of N known coarse grid node values (including required ghosts). Let G be the vector of S unknown fine grid node values. N <---------> X X X X X X ^ |g0| ^ X X X X X X | |f0| ^ |g1| | X X X X X X | |f1| | |g2| | X X X X X X | |f2| | |g3| | W = X X X X X X | S F= |f3| | N G = W*F = |g4| | S X X X X X X | |f4| | |g5| | X X X X X X | |f5| v |g6| | X X X X X X | |g7| | X X X X X X v |g8| v Will contain the same data type as intepolator.M if dtype is not passed, else W will be computed from user given dtype. Wr: np.ndarray Reduced W that exludes rightmost output points of ndimensional output vector. Pre computed weights to interpolate directly from coarse to inner fine grid. Let F be the vector of N known coarse grid node values (including required ghosts). Let G be the vector of GR unknown fine inner grid node values (see gr attribute). N <---------> X X X X X X ^ |g0| ^ X X X X X X | |f0| ^ |g1| | X X X X X X | |f1| | |g2| | X X X X X X | |f2| | |g3| | Wr = X X X X X X | GR F= |f3| | N G = W*F = |g4| | GR X X X X X X | |f4| | |g5| | X X X X X X | |f5| v |g6| | X X X X X X | |g7| | X X X X X X v |g8| v Same data type as W. """ check_instance(interpolator, PolynomialInterpolator) check_instance(grid_ratio, tuple, values=int, size=interpolator.dim, minval=1) p = interpolator.p n = interpolator.n P = interpolator.P N = interpolator.N M = interpolator.M ghosts = interpolator.ghosts gr = grid_ratio GR = np.prod(gr, dtype=np.int32) del grid_ratio s = tuple(gri + 1 for gri in gr) S = np.prod(s, dtype=np.int32) dim = interpolator.dim key = ("PolynomialSubgridInterpolator", interpolator.key, gr, str(dtype)) self.p = p self.P = p self.s = s self.S = S self.n = n self.N = N self.gr = gr self.GR = GR self.ghosts = interpolator.ghosts self.dim = dim self.interpolator = interpolator self.key = key cache_file = interpolator.cache_file() try: data = load_data_from_cache(cache_file, key) if data is not None: W = data if _check_matrix(W): self.W = data return except Exception as e: msg = f"Failed to load data from cache because:\n{e}" warnings.warn(msg, HysopCacheWarning) X = tuple( np.asarray(tuple(sm.Rational(j, gr) for j in range(0, si))) for i, (gr, si) in enumerate(zip(gr, s)) ) V = np.vander(X[0], N=p[0], increasing=True) for i in range(1, dim): Vi = np.vander(X[i], N=p[i], increasing=True) V = np.multiply.outer(V, Vi) even_axes = tuple(range(0, V.ndim, 2)) odd_axes = tuple(range(1, V.ndim, 2)) axes = even_axes + odd_axes V = np.transpose(V, axes=axes).copy() assert V.shape[:dim] == s assert V.shape[dim:] == p V = V.reshape((S, P)) W = V.dot(M) assert W.shape == (S, N) if dtype is not None: W = W.astype(dtype) update_cache(cache_file, key, W) self.W = W def __call__(self, F): return self.W.dot(F.ravel()).reshape(self.s)
[docs] def generate_subgrid_restrictor(self): return PolynomialSubgridRestrictor(subgrid_interpolator=self)
@property def Wr(self): assert self.W.shape == (self.S, self.N) view = (slice(0, -1, None),) * self.dim + (slice(None, None, None),) * self.dim Wr = self.W.reshape(self.s + self.n)[view].reshape(self.GR, self.N).copy() return Wr def __hash__(self): objs = (self.interpolator, self.gr) return hash(objs)
[docs] class PolynomialSubgridRestrictor: def __init__(self, subgrid_interpolator): """ Create a PolynomialSubgridRestrictor from a PolynomialSubgridInterpolator. Parameters ---------- subgrid_interpolator: PolynomialSubgridInterpolator Interpolant used to compute restrictor weights. Attributes ---------- g: tuple of int Corresponds to (n+1)*s G: int Corresponds to g[0]*g[1]*...g[dim-1] R: np.ndarray Restrictor weights. origin: tuple of int Origin of the generated stencil. Rr: np.ndarray Restrictor weights excluding leftmost and rightmost points. ghosts: Corresponds to origin - 1, which is also Rr origin. """ check_instance(subgrid_interpolator, PolynomialSubgridInterpolator) dim = subgrid_interpolator.dim n = subgrid_interpolator.n s = subgrid_interpolator.s gr = subgrid_interpolator.gr GR = subgrid_interpolator.GR W = subgrid_interpolator.W g = tuple(ni * gri + 1 for (ni, gri) in zip(n, gr)) G = np.prod(g, dtype=np.int64) assert all(gi % 2 == 1 for gi in g) origin = tuple(gi // 2 for gi in g) gvals, gvars = tensor_symbol("g", g, origin) I = 0 for idx in np.ndindex(*gr): mask = tuple( slice(i, max(2, i + ni * gri), gri) for (i, ni, gri) in zip(idx, n, gr) ) target = tuple(gri - i for (i, gri) in zip(idx, gr)) F = gvals[mask] Ii = W.dot(F.ravel()).reshape(s)[target] I += Ii R = np.ndarray(shape=g, dtype=object) for idx in np.ndindex(*g): R[idx] = I.coeff(gvals[idx]) view = (slice(1, -1, None),) * dim Rr = R[view] ghosts = tuple(oi - 1 for oi in origin) self.g = g self.G = G self.R = R self.Rr = Rr self.origin = origin self.ghosts = ghosts self.n = n self.s = s self.GR = GR self.gr = gr self.subgrid_interpolator = subgrid_interpolator
if __name__ == "__main__": np.set_printoptions( precision=4, linewidth=1e8, threshold=1e8, formatter={"float": lambda x: f"{x:+0.3f}"}, ) # 2D tests grid_ratio = (2, 2) F = [[1, 1], [1, 1]] F = np.asarray(F) print("Solving bilinear...") PI = PolynomialInterpolator(dim=2, deg=1, fd=2, verbose=False) GI0 = PI.generate_subgrid_interpolator(grid_ratio=grid_ratio) GI1 = PI.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("Bilinear (Rational)") print(GI0(F)) print() print("Bilinear (np.float64)") print(GI1(F)) print() F = [[0, 0, 0, 0], [0, 1, 1, 0], [0, 1, 1, 0], [0, 0, 0, 0]] F = np.asarray(F) print("Solving bicubic2...") PI0 = PolynomialInterpolator(dim=2, deg=3, fd=2, verbose=False) GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("Bicubic (FDC2)") print(GI0(F)) print() print("Solving biquintic2...") PI1 = PolynomialInterpolator(dim=2, deg=5, fd=2, verbose=False) GI1 = PI1.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("Biquintic (FDC2)") print(GI1(F)) print() F = [[0, 1, 1, 0], [0, 1, 1, 0]] F = np.asarray(F) print("Solving linear/cubic...") PI0 = PolynomialInterpolator(dim=2, deg=(1, 3), fd=2, verbose=False) GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("Linear/Cubic (FDC2)") print(GI0(F)) print() F = [[0, 0], [1, 1], [1, 1], [0, 0]] F = np.asarray(F) print("Solving cubic/linear...") PI0 = PolynomialInterpolator(dim=2, deg=(3, 1), fd=2, verbose=False) GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("Cubic/Linear (FDC2)") print(GI0(F)) print() F = [ [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0], [0, 0, 1, 1, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], ] F = np.asarray(F) print("Solving bicubic4...") PI0 = PolynomialInterpolator(dim=2, deg=3, fd=4, verbose=False) GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("Bicubic (FDC4)") print(GI0(F)) print() print("Solving biquintic4...") PI1 = PolynomialInterpolator(dim=2, deg=5, fd=4, verbose=False) GI1 = PI1.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("Biquintic (FDC4)") print(GI1(F)) print() print("Solving biseptic2...") PI2 = PolynomialInterpolator( dim=2, deg=7, fd=2, verbose=False, approximative=(not has_flint) ) GI2 = PI2.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("Biseptic (FDC2)") print(GI2(F)) print() print("Solving binonic2...") PI3 = PolynomialInterpolator( dim=2, deg=9, fd=2, verbose=False, approximative=(not has_flint) ) GI3 = PI3.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("Binonic (FDC2)") print(GI3(F)) print() print("Solving septic2/nonic2 ...") PI4 = PolynomialInterpolator( dim=2, deg=(7, 9), fd=2, verbose=False, approximative=(not has_flint) ) GI4 = PI4.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("septic/nonic (FDC2)") print(GI4(F)) print() print("Solving septic2/quintic4 ...") PI5 = PolynomialInterpolator( dim=2, deg=(7, 5), fd=(2, 4), verbose=False, approximative=(not has_flint) ) GI5 = PI5.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print("septic/nonic (FDC2/FDC4)") print(GI5(F)) print() # 3D test grid_ratio = (2, 2, 2) print("Solving trilinear...") PI = PolynomialInterpolator(dim=3, deg=1, fd=2, verbose=False) GI0 = PI.generate_subgrid_interpolator(grid_ratio=grid_ratio) print() print("Solving tricubic2...") PI0 = PolynomialInterpolator(dim=3, deg=3, fd=2, verbose=False) GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print() print("Solving triquintic2...") PI0 = PolynomialInterpolator( dim=3, deg=5, fd=2, verbose=False, approximative=(not has_flint) ) GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64) print()